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INTRODUCTION 


The  USAF  Phillips  Laboratory  Scholar  Program  is  sponsored  by  the  Phillips 

Laboratory  Geophysics  Directorate.  It  is  conducted  by  the  Northeast  Consortium  for 
Engineering  Education  (NCEE).  This  program  provides  research  opportunities  for 

qualified  engineers  and  scientists  who  have  received  their  PhD  degrees,  or 

equivalent,  from  technical  programs  at  U.S.  universities  or  technical  institutions. 

These  opportunities  consist  of  research  appointments  of  12  months  duration  with  the 

Geophysics  Directorate,  located  at  Hanscom  AFB  near  Boston,  Massachusetts,  or  the 
Advanced  Weapons  and  Survivability  Directorate  located  at  Kirtland  Air  Force  Base, 
New  Mexico. 


PROGRAM  DESCRIPTION 

The  Air  Force  Geophysics  Directorate  has  initiated  the  Phillips  Laboratory 
Scholar  Program  to  broaden  the  direct  participation  of  qualified  researchers  in 
Phillips  Laboratory  research  programs.  This  program  has  provided  research 

opportunities  for  selected  engineers  and  scientists  holding  a  doctoral  degree  to  work 
at  the  Air  Force  Phillips  Laboratory  for  a  12  month  research  period;  the  appointment 

may  be  extended  for  a  second  term. 

To  be  eligible,  candidates  must  have  a  PhD  or  equivalent  experience  in  an 
appropriate  technical  field.  Scholars  located  at  Hanscom  AFB  will  be  selected 
primarily  from  such  basic  and  applied  science  fields  as  physics,  particularly 
geophysics  and  atmospheric  physics,  chemistry  computer  science,  and  engineering. 

Scholars  located  at  Kirtland  AFB  will  be  involved  in  pulsed  power,  high-density  and 
high-energy  plasmas,  high-power  microwave  sources,  electromagnetic  applications 
and  effects,  satellite  survivability  and  vulnerability  analysis,  and  related 

experimental,  theoretical,  and  computational  method  topics.  Applicants  must  be  U.S. 

citizens  who  are  recent  recipients  of  the  PhD  (within  the  last  five  years). 

NCEE  supports  equal  opportunity/affirmative  action.  All  qualified  applicants 
will  receive  consideration  without  regard  to  race,  color,  religion,  sex,  or  national 

origin. 


PROGRAM  OBJECTIVES 

(1)  To  provide  a  productive  means  for  scientists  and  engineers  holding 
PhD  degrees  to  participate  in  research  at  the  Air  Force  Phillips 
Laboratory; 

(2)  To  stimulate  continuing  professional  association  among  the 
Scholars  and  their  professional  peers  in  the  Air  Force; 

(3)  To  further  the  research  objectives  of  the  United  States  Air  Force; 

(4)  To  enhance  the  research  productivity  and  capabilities  of  scientists 
and  engineers  especially  as  these  relate  to  Air  Force  technical 

interests. 

This  report  covers  the  period  1  July  1997  through  28  May  1998.  Four  scholars 
participated  in  the  Phillips  Laboratory  Scholar  Program  during  this  period  with  two 
completing  the  program. 
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Sean  J.  Carey 

My  tenure  as  a  Geophysics  Scholar  began  July  15,  1997  and  ends  April  21,  1998.  During  this  time, 
I  have  worked  on  two  separate  projects,  an  investigation  of  grain  destruction  around  mid  to  late  B 
stars  and  the  discovery  and  investigation  of  the  MSX  mid-infrared  dark  clouds.  This  final  report 
contains  three  sections.  The  first  two  sections  discuss  my  two  projects  with  a  brief  introduction,  a 
summary  of  important  results  and  a  short  description  of  the  work  completed.  The  final  section  of 
the  report  lists  the  papers  published,  presentations  given,  and  meetings  attended  during  my  tenure. 

Grain  Destruction  Around  B  Stars: 


Work  is  in  collaboration  with  Drs.  Frank  O.  Clark  and  Russell  F.  Shipman  of  VSBC/AFRL  and  Dr. 
Rob  Assendorp  of  the  Astrophysikalishches  Institut  Potsdam. 

The  exact  nature  of  the  composition  and  size  distribution  of  interstellar  dust  grains  is  a  current  and 
very  important  area  of  research.  Interstellar  dust  grains  are  10  A  to  1  pm  sized  solids  composed  of 
silicate  and  carbonaceous  compounds.  The  properties  of  interstellar  dust  grains  are  inferred  from 
the  visible  and  UV  starlight  they  absorb  and  scatter  (the  combination  of  absorption  and  scattering  is 
called  extinction)  and  re-radiate  in  the  infrared.  We  investigated  the  infrared  emission  and 
ultraviolet  extinction  of  interstellar  dust  in  which  main  sequence  B  stars  are  embedded  to  determine 
how  the  enhanced  radiation  field  near  the  star  affects  the  interstellar  dust.  In  particular,  we  studied 
the  properties  of  the  smallest  dust  grains  (<  0.1  pm),  as  they  are  the  most  affected  by  the  presence 
of  an  enhanced  radiation  field.  The  smallest  dust  grains  are  heated  stochastically  by  UV  photons 
and  emit  primarily  at  wavelengths  of  12-25  pm.  In  addition,  the  UV  extinction  is  more  sensitive  to 
variations  in  the  number  and  physical  properties  of  these  smaller  dust  grains. 

Our  study  was  motivated  by  previous  work  of  Shipman  (1995)  and  Shipman  &  Clark  (1995)  who 
studied  dust  emission  in  HII  regions.  They  determined  that  12  and  25  pm  (small  grain)  emission 
declined  relative  to  the  emission  at  100  pm  (due  primarily  to  larger  grains)  as  the  radiation  field 
increased.  They  interpret  the  relative  decline  in  12  and  25  pm  emission  to  the  destruction  of  small 
grains  by  the  high  flux  of  UV  photons  from  the  exciting  star  of  the  HU  region.  We  extend  this 
prior  work  by  studying  the  lower  UV  flux  regions  around  B  stars.  Specifically,  we  are  interested  in 
determining  the  hardness  of  the  radiation  field  necessary  to  produce  12  and  25  pm  emission  deficits 
and  using  the  results  to  determine  the  structure  of  the  small  grains.  Our  preliminary  investigation 
(Carey  et  al.  1997)  determined  that  the  12  and  25  pm  emission  was  deficient  in  the  region  around  B 
stars. 


Methodology: 

The  sample  set  of  stars  studied  was  taken  from  the  list  of  stars  embedded  in  infrared  cirrus  clouds 
complied  by  Gaustad  &  Van  Buren  (1993).  From  this  list,  we  selected  bright,  main  sequence  stars 
of  spectral  types  B1.5-B9  within  1000  parsecs.  The  IRAS  ISSA  images  (Wheelock  et  al.  1994)  at 
12,  25,  60  and  100  pm  were  examined  for  emission  around  each  star.  We  restricted  our  sample  set 
to  the  20  stars  with  the  least  confused  infrared  backgrounds.  The  infrared  emission  was  averaged 
in  annuli  around  each  star  to  construct  an  emission  profile  as  a  function  of  angular  distance  from 


the  star  at  12,  25,  60  and  100  pm.  A  background  was  subtracted  from  each  annulus  to  account  for 
foreground  emission  and  the  uncertainty  in  the  DC  offset  level  of  the  ISSA  images.  Detector 
hysteresis  and  AC/DC  responsivity  variations  were  investigated  and  determined  not  to  affect  the 
quality  of  the  data.  Color  ratios  were  calculated  from  the  averaged  emission  profiles.  The  color 
ratio  profiles  were  compared  with  the  grain  emission  models  of  Desert,  Boulanger  and  Puget  (DBP; 
1990).  The  model  results  were  converted  to  angular  distance  from  the  star  and  ratioed  to  produce 
color  profiles  for  comparison. 

We  measured  the  UV  extinction  toward  four  of  the  stars  using  UV  photometry  from  the 
Astronomical  Netherlands  Satellite.  We  used  the  UV  color  excesses  published  by  Savage  et  al. 
(1985).  We  fit  the  observed  extinction  curves  with  the  model  of  Clayton,  Cardelli  and  Mathis  to 
derive  the  ratio  of  selective-to-total  extinction  (Rv)  for  the  lines  of  sight  towards  the  stars  studied. 

Results: 

For  all  sources,  the  infrared  emission  at  12, 25, 60  and  100  pm  increases  as  the  star  is  approached. 
The  60/100  color  is  well  modeled  by  DBP  rising  from  a  uniform  background  value  of  0.2  to  a  peak 
at  the  position  of  the  star.  Surprisingly,  the  rise  in  25  pm  emission  begins  at  a  larger  distance  from 
the  star  than  the  60  pm  increase.  It  appears  that  the  25  pm  emitters  are  more  responsive  than  the 
rest  of  the  grains  to  slight  changes  in  the  radiation  field.  Five  stars  have  declining  12/100  and 
25/100  colors  as  the  star  is  approached,  which  is  in  direct  conflict  with  the  DBP  model.  The 
25/100  decline  is  more  pronounced  than  the  12/100  decline.  The  12/100  and  25/100  color  deficits 
are  associated  with  far-UV  extinctions  lower  than  the  standard  interstellar  value.  Destruction  of 
small  grains  in  the  increased  radiation  field  near  the  star  is  necessary  to  explain  the  correspondence 
of  12/100  and  25/100  color  deficits  with  low  far-UV  extinction.  The  particulars  of  the  color 
variation  imply  that  the  25  pm  emitters  are  preferentially  destroyed.  Our  investigation  suggests 
that  the  carriers  of  the  25  pm  emission  are  loosely  bonded  aggregates  of  carbonaceous  compounds 
such  as  polycyclic  aromatic  hydrocarbons.  These  aggregates  can  be  photodissociated  by  the 
moderate  enhancement  in  the  ultraviolet  radiation  field  found  around  main  sequence  B  stars.  Our 
analysis  is  consistent  with  the  conglomerate  grain  model  of  Clark  et  al.  (1995). 

References: 

Clayton,  G.  C.,  Cardelli,  J.  A.,  &  Mathis,  J.  S.  1988,  ApJ,  329,  L33 
Carey,  S.  J.,  Shipman,  R.  F.,  Clark,  F.  O.,  &  Assendorp,  R.  1997,  ApJ,  479,  303 
Clark,  F.  O.,  Shipman,  R.  F.,  Assendorp,  R.,  Kester,  D.,  &  Egan,  M.  P.  1995,  Plan.  Sp.  Sci.,  43, 
1353 

Desert,  F.-X.,  Boulanger,  F.,  &  Puget  J.  L.  1990,  A&A,  237,  215 
Gaustad,  J.  E.,  &  van  Buren,  D.  1993,  PASP,  105, 1127 
Savage,  B.  D.,  Massa,  D.,  Meade,  M.,  Wesselius,  P.  1985,  ApJS,  59, 397 
Shipman,  R.  F.  1994,  PhD.  Thesis,  University  of  Wyoming 
Shipman,  R.  F.,  &  Clark,  F.  O.  1994,  ApJ,  422, 153 

Wheelock,  S.  L.  et  al.  1994,  IRAS  Sky  Survey  Atlas  Explanatory  Supplement  (JPL  Publ.  94-11; 
Pasadena:  JPL) 
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MSX  Infrared-Dark  Clouds: 


Work  is  in  collaboration  with  Drs.  Michael  P.  Egan,  Frank  O.  Clark,  Stephan  D.  Price,  Russell  F. 
Shipman  of  VSBC/AFRL  and  Thomas  A.  Kuchar  of  Dynamic  Research  Corporation 

The  SPIRIT  III  infrared  telescope  aboard  the  Midcourse  Space  Experiment  (MSX)  observed  the 
Galactic  plane  with  high  spatial  resolution  (18")  and  sensitivity  in  five  wavelength  bands  between  4 
and  25  |xm  (Price  1995).  Quite  surprisingly,  the  bright  mid-infrared  emission  of  the  Galactic  plane 
is  punctuated  by  dark  patches.  These  dark  regions  are  objects  foreground  to  the  majority  of  the 
emission  with  high  mid-infrared  opacities.  Motivated  by  the  discovery  of  approximately  2000  of 
these  clouds  in  a  cursory  inspection  of  initially  processed  images  of  the  Galactic  plane,  we  have 
initiated  a  multi- wavelength  study  of  these  infrared-dark  regions  (hereafter  called  MSX  dark 
clouds)  to  better  understand  their  physical  properties. 

Methodology: 

One  of  the  goals  of  the  MSX  Galactic  plane  survey  is  to  produce  well  calibrated  images  of  the 
Galactic  plane  in  wavelength  bands:  6.8-10.8  (A),  4.22-4.36  (Bl),  4.24-4.45  (B2),  11.1-13.2  (C), 
13.5-15.9  (D)  and  18.2-25.1  (E)  pm.  To  achieve  this  result,  a  number  of  data  artifacts  and 
calibration  issues  must  be  handled.  Accurate  calibration  of  the  MSX  data  is  needed  to  calculate  the 
mid-infrared  extinction  across  the  MSX  dark  clouds.  A  crucial  task  is  the  proper  subtraction  of  the 
instrumental  dark  current.  The  initial  dark  current  correction  provided  is  not  sufficient  for 
astronomical  purposes.  I  have  been  involved  in  attempts  to  produce  better  corrections  by  modeling 
the  measured  dark  current  measurements  in  more  detail.  This  work  is  ongoing,  but  initial  results 
show  that  it  is  feasible  to  more  accurately  model  the  temporal  variations  of  the  dark  current. 

In  addition  to  a  more  detailed  study  of  the  mid-infrared  data,  the  MSX  dark  clouds  have  been 
observed  in  millimeter- wave  rotational  transitions  of  H^CO.  H2CO  is  a  tracer  of  the  dense 
molecular  gas  that  should  be  associated  with  the  large  column  of  dust  producing  the  observed  mid- 
IR  extinction.  We  observed  ten  clouds  using  the  NRAO  12  Meter  telescope  at  Kitt  Peak,  Arizona. 
Where  possible,  we  have  mapped  the  MSX  dark  clouds  in  the  2 12- 111  lfr>e  of  H2CO.  By  observing 
several  different  transitions  of  ortho  and  para-H2CO  and  the  isotopic  substitutions,  H213CO  and 
HDCO,  we  have  probed  the  gas  density,  kinetic  temperature,  column  density  and  H2CO  abundance 
of  these  clouds.  Line  properties  have  been  determined  by  the  fitting  of  Gaussian  profiles  and/or 
direct  calculation  of  line  moments.  The  observed  transitions  have  been  modeled  using  a  large 
velocity  gradient  (LVG)  statistical  equilibrium  code  (Mangum  &  Wootten  1993).  Distances  to 
some  of  the  clouds  have  been  determined  from  the  kinematic  information  provided  by  the  observed 
lines  and  a  model  of  the  Galactic  rotation  curve  (Clemens  1985). 

The  far-infrared  emission  of  the  MSX  dark  clouds  has  been  studied  using  the  IRAS  ISSA  plates 
(Wheelock  et  al.  1994),  IRAS  Galaxy  Atlas  (Cao  et  al.  1997)  and  IRAS  Groningen  GEISHA  high 
resolution  images  (Assendorp  et  al.  1995).  We  could  not  attribute  far-infrared  emission  above  the 
local  background  (1000  MJy/sr)  to  any  of  the  clouds  observed  in  H2CO.  An  upper  limit  to  the  dust 
temperature  for  the  clouds  was  set  by  modeling  the  expected  emission  using  grain  absorption 
coefficients  from  Mathis  &  Whiffen  (1989)  and  assuming  that  the  emission  was  due  to  an  optical 
thin  column  of  0.1  Jim  sized  grains.  We  have  proposed  submillimeter  continuum  observations 
using  the  SCUBA  camera  on  the  James  Clerk  Maxwell  Telescope  to  search  the  clouds  for 
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embedded  protostellar  objects  and  to  defmitively  determine  the  column  density  and  temperature  of 
the  interstellar  dust  in  these  clouds. 

Results: 

All  ten  clouds  observed  in  H2CO  2)2-1 11  were  detected,  confirming  that  the  MSX  dark  clouds 
contain  dense  molecular  gas.  The  morphologies  of  the  spectral  line  emission  vary  from  globule 
shaped  to  filamentary  and  closely  follow  the  mid-infrared  extinction.  We  determined  distances 
between  1.0  to  8.5  kiloparsecs  to  the  clouds.  Also,  the  clouds  are  kinematically  associated  with  HII 
regions  (Kuchar  &  Clark  1997)  suggesting  that  they  are  associated  with  star  forming  regions.  The 
linear  diameters  of  the  MSX  dark  clouds  are  0.2-10  parsecs.  The  molecular  column  densities  are 
between  1023  and  1025  cm2  with  estimates  from  mid-infrared  extinction  and  H2CO  emission 
agreeing. 

From  LVG  modeling,  we  estimate  that  the  gas  densities  of  the  clouds  are  of  order  106  cm’3  and  the 
gas  kinetic  temperatures  are  10-15  K.  The  upper  limit  to  the  dust  temperature  from  the  lack  of  far- 
infrared  emission  is  25  K.  The  gas  and  dust  temperatures  should  be  equilibrated  at  the  densities  of 
the  MSX  dark  clouds.  The  abundance  of  H2CO  relative  to  H2  is  —1.0  x  10'10  and  is  significantly 
depleted  compared  to  more  diffuse  molecular  clouds.  It  is  likely  that  H2CO  is  depleted  by 
accretion  of  gas-phase  metals  onto  the  dust  grains  in  the  cold,  dense  MSX  dark  clouds. 

The  MSX  dark  clouds  have  nearby  star-forming  regions  but  do  not  have  any  known  star-forming 
tracers  (HE  regions,  masers  or  IR  point  sources)  within  their  boundaries.  Given  the  high  densities 
of  these  objects,  it  is  unlikely  that  they  are  stable  against  collapse  and  subsequent  star  formation 
despite  the  lack  of  star  formation  tracers.  Indeed,  the  observed  molecular  line  profiles  have 
pronounced  wings,  which  suggest  some  sort  of  internal  motions  in  these  clouds.  It  is  likely  that  the 
MSX  dark  clouds  contain  internal  structure  in  a  pre-protostellar  state  and  may  contain  class  0 
protostars.  Further  observations  (particularly  in  the  submillimeter)  will  be  able  to  more  accurately 
characterize  the  star-forming  properties  of  these  objects. 

The  MSX  Galactic  plane  survey  has  identified  a  new  population  of  objects.  These  objects  are 
extremely  large,  cold,  dense  molecular  cores  associated  with  star  forming  regions,  but  clearly  in  a 
protostellar  or  pre-protostellar  state  themselves.  Further  study  of  these  exciting  objects  should 
provide  crucial  information  to  understanding  the  earliest  stages  of  the  star  formation  process. 
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Physical  Properties  of  the  MSX  Galactic  Infrared-Dark  Clouds",  submitted  to  the 
Astrophysical  Journal 

Carey,  Sean  J.,  Shipman,  R.  F.,  &  Clark,  F.  O.,  "Deficits  of  25  pm  Carriers  Around  B  Stars",  in 
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Presentation  given: 
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"Infrared  Dark  Clouds  in  the  Milky  Way",  poster  presented  at  the  191st  meeting  of  the 
American  Astronomical  Society,  Washington,  DC 

Proposals  submitted: 

"Molecular  Search  for  MSX  Dark  Clouds",  proposal  for  NRAO  12  meter  observing  time, 
observations  conducted  November  1-6,  1997 

NRAO  12  meter  unassigned  time  request  for  January  1998,  observations  conducted  January  12, 
1998 

NRAO  12  meter  unassigned  time  request  for  March  1998,  observations  conducted  March  27,  1998 

"Study  of  Physical  Properties  of  MSX  Dark  Clouds",  proposal  for  NRAO  12  meter  observing  time, 
awarded  time  May  16-20,  1998 

"SCUBA  survey  of  MSX  Infrared-Dark  Clouds",  proposal  for  James  Clerk  Maxwell  Telescope 
observing  time,  submitted 

"Investigation  of  Dust  Chemistry  in  Cold  Molecular  Cores",  proposal  for  NRAO  12  meter 
observing  time,  submitted 


Travel: 

Observed  at  the  NRAO  12  meter  telescope,  Kitt  Peak,  Arizona,  November  1-6, 1997 

Attended  191st  meeting  of  the  American  Astronomical  Society,  Washington,  DC,  January  6-10, 
1998 
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INTRODUCTION 

I  describe  in  this  document  the  summary  of  accomplishments  during 
my  two-year  tenure  as  a  NCE£  Geophysics  Scholar  contracting  at  Air 
Force  Research  Laboratories  (formerly  Phillips  Laboratory).  Hanscom 
AFB,  Massachusetts.  I  describe  the  work  done  which  led  to  submission 
of  papers  for  refereed  publication,  and  I  outline  the  results  of  my  at¬ 
tendance  at  professional  conferences  during  my  tenure  as  a  Geophysics 
Scholar. 

WORK  SUBMITTED  FOR  PUBLICATION 

1.  I  researched  the  kinematics  of  small,  isolated,  structurally  simple 
star-forming  clouds  in  our  Galaxy  called  ‘:Bok  globules.’’  Fifteen  small 
Bok  globules  showing  no  signs  of  current  star  formation  were  mapped 
at  high  spatial  and  especially  spectral  resolution  in  the  (.7=1—7 
0)  rotational  lines  of  CO,  uCO,  and  C180,  using  the  14  meter  radio 
telescope  of  the  Five  College  Radio  Astronomy  Observatory.  Maps  were 
made  with  the  fifteen-element  3  mm  array  receiver  QUARRY  and  the 
15  x  1024  channel  autocorrelator  spectrometer  FAAS.  From  120  to  360 
positions  per  globule,  sampled  with  half-beam  spacing,  were  observed 
in  the  13CO  line,  while  30  to  60  full-beam-sampled  Cl80  and  CO 
positions  per  globule  were  observed.  l3CO  and  Cl80  were  observed 
with  a  velocity  resolution  ofless  than  0.007  kms-1  channel -1;  CO  was 
observed  with  0.013  kms"1  channel-1  resolution. 

Gaussian  fitting  of  the  emission  lines  was  used  to  establish  mean 
radial  velocities  and  uncertainties.  Each  globule  radial  velocity  dis¬ 
tribution  on  the  sky  was  fit  to  a  plane  (solid  body  rotation)  to  yield 
mean  velocity  gradients  with  position,  and  rotation  axis  directions. 
The  globules  were  found  to  be  rotating  at  rates  more  than  an  order 
of  magnitude  faster  than  velocity  shifts  attributable  to  local  differ¬ 
ential  Galactic  rotation.  For  globule  assumed  mean  distances  of  600 
pc.  the  gradients  range  from  0.089 kms-1  pc-1  (a/~-3  x  10“ 15  s-1)  to 
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0.950  kms”1  pc”1  (a>~ 3  x  10”14  s”1).  Half  show  gradients  less  than 
about  0.3 kms”1  pc-1,  and  half  show  gradients  greater  than  about 
0.6 kms”1  pc”1,  distributed  in  a  distinctly  bimodal  fashion. 

Detailed  examination  of  the  globule  rotation  curves  indicated  that 
the  kinematics  of  ten  of  fifteen  globule  cores  are  well-approximated  by 
solid-body  rotation.  Differential  rotation  and  shearing  motions  due  to 
external  influences  (ram  pressure  stripping  and/ or  bcrw  shocks)  are  also 
seen. 


2.  I  researched  the  stellar  content  of  some  of  the  Bok  Globules 
described  in  item  1  (and  of  others).  Using  the  Infrared  Space  Obser¬ 
vatory  (ISO),  we  took  ISOCAM  images  at  6  microns  and  15  microns 
wavelength,  in  order  to  confirm  previous  implications  from  Infrared  As¬ 
tronomy  Satellite  (IRAS)  data  that  said  globules  were  star-free.  ISO 
images  provided  much  more  sensitivity  noth  which  to  make  that  dec¬ 
laration.  Of  the  dozen  or  so  globules  studied  for  this  paper,  all  but 
one  were  decided  to  be  definitively  star-free.  The  remaining  cloud  will 
need  further  study  to  determine  its  characteristics. 

3.  I  researched  the  magnetic  field  configurations  around  the  same 
Bok  globules  described  in  item  1.  Fifteen  small  Bok  globules  were 
probed  using  a  CCD  imaging  polarimeter  in  order  to  create  detailed 
maps  of  the  magnetic,  fields  associated  with  these  clouds.  Stars  as 
faint  as  20th  mag  at  F-band  were  measured  polarimetrically  with  un¬ 
certainties  less  than  1%.  Sky  transmission  variations  were  minimized 
via  a  system  of  synchronous  polaroid  rotation  and  bidirectional  charge 
shifting.  In  all,  more  than  1000  stars  behind  the  periphery  of  these 
globules  were  accurately  analyzed  polarimetrically.  The  large-scale  (1 
-  2  pc)  magnetic  field  patterns  around  the  SBGs  may  be  classified  into 
four  types:  Type  I.  two  narrowly  dispersed  patterns  which  overlap  spa¬ 
tially;  Type  II,  one  widely  dispersed  pattern  and  one  narrowly  dispersed 
pattern  which  overlap  spatially;  Type  III,  primarily  uniform;  and  Type 
TV,  enmetary  (polarization  vectors  appear  to  wrap  around  the  globule). 
Multiple  Gaussian  fitting  of  the  polarization  position  angle  histograms 
results  in  la  dispersions  which  average  around  10°,  indicating  one  or 
more  predominantly  uniform  field  directions  along  the  line  of  sight. 
The  widely  dispersed  primary  peaks  of  patterns  of  Type  II  and  IV  are 
likely  due  to  changes  in  the  line-of-sight  field  direction,  as  evidenced 


by  the  smooth,  changes  in  position  angles  with  polarization.  The  sec¬ 
ondary  peaks  of  patterns  of  Types  I,  II,  and  III  are  typically  produced 
by  highly  polarized,  faint  background  stars;  likely  the  secondary  peaks 
trace  magnetic  field  conditions  which  are  not  local  to  the  SBGs. 

Among  the  radial  distributions  of  SBG  mean  polarization  and  po¬ 
sition  angle,  there  is  no  overall  trend.  However,  SBGs  with  patterns  of 
Type  II  or  III  tend  to  have  stars  with  higher  mean  polarizations  nearer 
to  their  centers,  usually  within  V  to  4'.  Among  the  corresponding 
distributions  of  eight  45u  azimuth  sectors  around  each  SBG,  most  of 
the  distributions  indicate  two  or  more  distinct  zones  of  polarizations 
and  position  angles.  Such  zones  are  most  distinct  around  SBGs  with 
field  patterns  of  Types  II  and  IIL  The  different  zones  likely  represent 
regions  with  different  line-of-sight  field  directions,  or  regions  with  mag¬ 
netic  fields  tangled  to  different  extents.  This  phenomenon  is  best  seen 
in  plots  of  mean  polarization  versus  dispersions  of  position  angles  in  the 
azimuthal  sectors  around  each  globule.  Almost  always,  regions  of  high 
mean  polarization  correspond  to  regions  of  low  dispersion  in  position 
angle,  and  vice  versa. 

Finally,  by  plotting  stellar  separations  against  differences  of  polar¬ 
ization  position  angles,  the  typical  SBG  was  found  to  have  a  magnetic 
field  structure  function  with  local  maxima  at  separations  of  several  ar- 
cminutes  (several  tenths  of  a  parsec:  the  mean  diameter  of  SBGs  in 
the  sample)  and  of  several  tenths  of  an  arcminute  (several  hundredths 
of  a  parsec:  a  possible  preferred  clump  size  in  SBG  cores). 

4.  I  helped  produce  a  high-resolution  catalog  of  our  Galactic  plane, 
at  wavelengths  of  12  and  25  microns.  We  used  a  maximum  entropy 
method  to  improve  resolution  on  Infrared  Astronomy  Satellite  (IRAS) 
data,  and  Cramer-Rao  bound  techniques  to  determine  noise  charac¬ 
teristics  and  accurate  fluxes.  Previous  catalogs  using  IRAS  data  have 
avoided  the  Galactic  plane  (from  ten  degrees  declination  on  either  side, 
and  out  to  100  degrees  from  the  Galactic  center),  because  the  many 
sources  in  this  region  led  to  confusion  effects  and  indeterminate  fluxes. 
This  work,  combined  with  the  similar  work  of  other  groups  at  60  and 
100  micron  wavelengths,  will  provide  valuable  spectral  information  tor 
a  variety  of  Galactic  and  extragalactic  sources,  ranging  from  stars  to 
nebulae.  Over  200,000  high-reliability  sources  were  catalogued  in  one 
wavelength  band  or  the  other,  or  in  both  bands.  Reliability  was  good 


down  below  0.2  Janskvs  in  flux. 


5.  I  researched  the  density  structures  of  the  same  Bok  globules 
described  in  items  1  and  3.  Fifteen  small  Bok  Globules  were  observed 
in  the  ( J  =  1  — *>  0)  rotational  lines  of  CO,  l3CO,  and  CldO,  using  the  14 
meter  radio  telescope  of  the  Five  College  Radio  Astronomy  Observatory 
(FCRAO).  Maps  were  made  with  the  fifteen-element  QUARRY  3  mm 
array  receiver  and  the  1024  channel  FAAS  autocorrelator  spectrometer. 
A  total  of  120  to  360  positions  per  globule,  sampled  with  half-beam 
spacing,  was  observed  in  the  l3CO  line,  mostly  over  2.5  by  3  arcmin 
grids,  while  30  to  60  full-beam-sampled  ClS0  and  CO  positions  per 
globule  were  observed.  13CO  and  Cl80  were  observed  with  a  velocity 
resolution  of  less  than  0.007  kms-1  channel-1;  CO  was  observed  with 
0.013  kms-1  channel-1  resolution. 

The  median  SBG  in  the  sample  has  an  apparent  volume  density 
profile  which  falls  off  as  ~  r-2*6,  significantly  steeper  than  found  in 
Yun  &  Clemens’  1991  sample  (of  which  the  majority  contained  YSO 
candidates)  whose  mean  dust  density  profile  fell  off  as  a  much  shallower 
~  r-1-6. 

Overall,  assuming  a  distance  of  600  pc.  the  sample  of  SBGs  is  char¬ 
acterized  as  containing  cores  of  size  0.3  pc.  SBGs  contain  about  10 
of  gas,  have  mean  H2  densities  of  a  few  xlO3  cm-3,  and  kinetic  temper¬ 
atures  of  ~  10  K.  Most  SBGs  are  near  virial  equilibrium.  Using  energy 
balance  arguments,  half  of  the  SBGs  may  be  quasi-statically  contract¬ 
ing,  and  all  of  these  globules  still  possess  envelopes.  The  strongest  as¬ 
sociation  with  cloud  contraction  is  the  existence  of  a  relatively  large  en¬ 
velope.  Rotation  is  in  most  cases  the  least  significant  source  of  support 
against  gravitational  contraction;  in  the  majority  of  cases  kinetic  en¬ 
ergy  owing  to  turbulent  motions  is  providing  the  most  support  against 
contraction. 

6.  I  initiated  PDSL  work  to  develop  a  web  site  devoted  to  astronomy 
and  space  science  and  technology  education. 

Summary:  Please  see  the  URL  http://spdcc.com/  kane/resume.html 


CONFERENCES  ATTENDED 


1.  My  tenure  at  the  Air  Force  Research  Laboratories  commenced 
with  the  publication  of  a  paper  given  at  the  ”Polarimetry  of  the  In¬ 
terstellar  Medium”  conference  in  Troy,  NY,  shortly  before  I  began  in 
my  position.  This  conference  was  devoted  both  to  the  advancement  of 
polarization  theories  and  to  the  dissemination  of  the  latest  results  in 
polarization  studies  of  the  space  between  stars  which  is  known  to  form 
new  stars.  My  paper  was  one  of  the  latter  observational  studies.  Re¬ 
visions  ensued,  leading  to  the  appearance  of  my  paper  “Investigating 
Correlations  in  the  Kinematics  and  Magnetic  Fields  of  Quiescent  Bok 
Globules”  in  ASP  Conference  Series  Vol.  97,  ed.  W.  G.  Roberge  &:  D. 
C.  B.  Whittet  (San  Francisco:  ASP),  p.  269  (1996). 

2.  In  October.  1996,  I  attended  the  Seventh  Astrophysics  Confer¬ 
ence  at  the  University  of  Maryland,  College  Park.  This  conference  was 
conceived  to  invite  discussion  on  the  current  state  of  star  formation 
theories  and  the  data  currently  available  to  confirm  or  rule  out  such 
theories.  As  such,  it  was  a  more  general  conference  than  the  Polarime- 
try  conference.  I  gave  the  paper  “Rotation  of  Starless  Bok  Globules,” 
as  a  prelude  to  my  published  work.  This  paper  can  be  found  in  A1P 
Conference  Series.  Vol.  393,  ed.  S.  H.  Holt  &  L.  G.  Mundy  (Woodbury, 
NY:  AIP),  p.  137  (1997). 

3.  In  January,  1997,  1  attended  the  189th  meeting  of  the  American 
Astronomical  Society,  held  in  Toronto,  Ontario,  Canada.  I  presented 
at  work  entitled  “The  Kinematics,  Stuctures,  and  Magnetic  Natures 
of  Starless  Bok  Globules” ,  which  summarized  my  Ph.  D.  dissertation 
and  the  continuing  work  I  had  done  as  a  Geophysics  scholar  to  improve 
upon  it. 

4.  In  January',  1998, 1  will  attend  the  191st  meeting  of  the  American 
Astronomical  Society,  held  in  Washington,  DC.  I  will  present  a  paper 
entitled  “Density  Structures  of  Starless  Bok  Globules,”  which  will  be 
a  detailed  look  at  research  commenced  for  my  Ph.  D.  dissertation  and 
■which  was  the  last  project  I  had  worked  on  as  a  Geophysics  Scholar. 
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During  the  aforementioned  contract  period,  I  have  utilized  a  high  temperature 
flowing  afterglow  (HTFA)  to  measure  rate  constants  for  ion-molecules  reactions  at  high 
temperatures  under  the  direction  of  Dr.  Albert  A.  Viggiano  at  the  Air  Force  Research 
Laboratory  in  the  Space  Vehicles  Directorate  (VSBP)  at  Hanscom  AFB,  MA.  Plasma 
chemistry  occurs  at  very  high  temperatures  and  has  a  critical  impact  on  Air  Force 
systems.  A  plasma  develops  around  re-entry  vehicles  that  can  interfere  with 
communications.  In  addition,  the  chemistry  of  the  ionosphere  occurs  at  temperatures 
upwards  of  2000  K.  Consequently,  accurate  modeling  of  ionospheric  processes 

requires  measurements  of  reactions  at  these  high  temperatures.  1  Also,  interest  in 
designing  enhanced  hypersonic  vehicle  combustors  that  employ  ion  chemistry  also 
necessitates  high  temperature  data  to  develop  similar  models. 2.3  Thus,  the  previously 

mentioned  HTFA  apparatus  has  been  designed  to  address  these  problems.^ 

The  high  temperature  flowing  afterglow  (HTFA)  operates  in  the  following 
manner.^  A  heated  filament  emits  energetic  electrons  that  transfer  energy  to  a  fast  flow 
of  helium  buffer  gas  through  collisions  that  produce  electronically  excited  helium  atoms. 
These  excited  atoms  transfer  their  energy  to  a  source  gas  introduced  downstream  from 
the  filament,  whereby  the  source  gas  becomes  ionized.  After  tens  of  centimeters,  the 
ions  equilibrate  to  the  flow  tube  temperature,  which  is  established  by  heating  the  tube 
with  a  commercially  available  furnace.  A  known  concentration  of  reactant  gas  is 
introduced  and  the  decrease  in  the  source  ion  signal  is  measured.  Plotting  the  decline 
of  the  ion  signal  as  a  function  of  reactant  concentration  over  a  known  reaction  time 
allows  the  rate  constant  to  be  calculated. 

Two  reactions  that  have  been  studied  are  the  charge  transfer  reactions  of  Ar+ 
with  02  and  CO  to  produce  02+  and  CO+,  respectively,  as  given  in  Equations  (1)  and  (2). 

Ar+|2P3/2j  +  02  ->  02  +  Ar  (1) 


1 


(2) 


Ar+|2P3/2)  +  CO  -»  CO+  +  Ar 

The  rate  constants  for  both  reactions  have  been  previously  measured  over  a  wide  range 
of  translational  energies  in  a  flow  drift  tube  at  room  temperature  by  Dotan  and 
Lindinger.5  The  effect  of  rotational  excitation  on  these  reactions  has  been  studied 
previously  in  our  lab  using  a  variable  temperature-selected  ion  flow  drift  tube  (VT- 
S1FDT).  However,  the  results  are  inconclusive.®  The  conditions  in  the  HTFA  produce 
substantial  rotational  as  well  as  vibrational  excitation  at  the  highest  temperatures. 
Therefore,  comparing  the  rate  constants  measured  in  the  HTFA  with  the  previous  drift 

tube  work  of  Dotan  and  Lindinger5  demonstrates  the  role  of  internal  energy  in  the 

charge  transfer  between  Ar+  with  O2  and  CO.^ 

Figure  1  shows  the  rate  constants  for  reactions  (1)  and  (2)  measured  in  the 
HTFA  plotted  against  the  sum  of  the  average  translational  and  rotational  energies. 

When  plotted  as  such,  the  HTFA  results  taken  at  temperatures  up  to  900  K  agree  well 
with  the  flow  drift  tube  results  of  Dotan  and  Lindinger5  given  by  the  open  symbols.  This 
agreement  indicates  that  rotational  and  translational  energy  is  equivalent  for  these 
reactions.  Below  900  K,  the  rate  constants  decrease  with  increasing  temperature.  The 
reaction  occurs  through  a  collision  complex  that  exists  for  a  few  rotational  periods  such 
that  the  energy  can  be  redistributed.  If  the  reaction  occurs  via  multiple  passes  through 
an  avoided  crossing  on  the  potential  energy  surfaces  and  the  lifetime  of  the  complex 
decreases  with  increasing  temperature,  the  complex  has  less  time  to  access  the 
crossing  seam.  The  transition  probability  for  crossing  the  seam  decreases,  in  turn 
decreasing  the  reaction  probability.  However,  the  rate  constants  increase  with 
increasing  temperature  above  900  K.  The  rate  constants  for  reaction  from  the 
vibrationally  excited  states  of  the  02  and  CO  reactants  can  be  calculated,  if  the  total  rate 
constant  observed  equals  the  sum  of  the  rate  constants  from  the  various  vibrational 
levels  weighted  by  their  statistical  populations.  We  also  assume  that  the  v  =0  and  1 

levels  both  react  at  the  drift  tube  rate  of  Dotan  and  Lindinger5  and  that  all  vibrational 
levels  above  v"=1  react  at  the  same  rate.  The  resulting  values  given  in  Figure  1  are 
independent  of  temperature  and  are  identical  to  the  collision  rate  constant.  This 
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observation  lends  supports  the  argument  that  the  v"=2  levels  are  responsible  for  the 
enhancement  in  the  rate  constants  at  high  temperatures,  which  has  been  observed  in 
other  systems  as  well.1  Populating  this  level  at  higher  temperatures  makes  it 
energetically  allowable  under  the  experimental  conditions  to  produce  charge  transfer 
products  in  an  excited  electronic  state,  a  process  that  is  quite  rapid.  Consequently,  the 

rate  constants  increase  when  this  reaction  path  becomes  available.5^ 

A  similar  reaction  studied  in  the  HTFA  is  the  charge  transfer  of  Ar+  with  C02 

shown  in  Equation  3  that  has  also  been  studied  in  a  flow  drift  tube.5 

Ar+|2P3/2j  +  C02  ->  C02  +  Ar  (3) 

The  results  for  reaction  (3)  are  shown  in  Figure  2  along  with  the  previous  flow  drift  tube 
results5  plotted  versus  the  average  total  energy,  i.  e„  the  sum  of  the  average  rotational, 
translational  and  vibrational  energies.  When  plotted  this  way,  the  results  from  the  two 
different  methods  agree  well  except  at  the  highest  few  temperatures.  The  reactivity 
consequently  depends  on  the  total  energy  available.  The  rate  constants  decrease  with 
increasing  temperature,  again  consistent  with  having  a  collision  complex  that  lives  for  a 
few  rotational  periods  such  that  the  energy  can  be  redistributed.  However,  at  the 
highest  temperatures,  the  rate  constants  obtained  in  the  HTFA  are  lower  than  those 
observed  in  the  flow  drift  tube.  If  the  lifetime  of  the  complex  decreases  with  increasing 
temperature,  the  complex  may  not  survive  for  enough  rotational  periods  at  the  highest 
temperatures,  thus  decreasing  the  reaction  rate.  In  this  case,  vibrational  excitation 
hinders  the  reaction  rate,  as  opposed  to  reactions  (1)  and  (2)  where  the  excitation 

enhances  it. 

A  third  reaction  examined  using  the  HTFA  is  the  charge  transfer  reaction 
between  02+  and  NO  which  is  important  in  the  chemistry  of  the  ionosphere.  The 
reaction,  given  by  Equation  (4)  is  important  in  maintaining  the  balance  of  ionization  in 
the  D  and  E  regions  from  70  to  150  km  in  altitude. 

oj  (x  2n  g)  +  no  -»  no+  +  o2  (4) 

Because  of  its  role  in  atmospheric  chemistry,  the  reaction  has  been  studied  as  a 
function  of  kinetic  energy  using  a  flow  drift  tube5  and  as  a  function  of  temperature  using 
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a  flowing  afterglow9  by  Lindinger,  et.  al.  The  flowing  afterglow  study,  however,  only 
includes  data  up  to  900  K.  The  current  experiments  extend  this  study  to  1400  K  using 
the  HTFA  apparatus.  Figure  3  shows  the  rate  constants  obtained  in  the  HTFA  plotted 
versus  temperature,  along  with  the  earlier  flow  drift  tube[Lindinger,  1974  #478  and 


flowing  afterglow9  measurements.  All  of  the  experiments  show  that  the  rate  constants 
are  independent  of  temperature,  giving  a  value  that  is  about  half  of  the  collision  rate 
constant.  The  current  results  may  show  a  slight  increase  with  increasing  temperature  at 
the  highest  temperatures.  Considering  the  15%  error  in  the  HTFA  measurements,  the 
rate  constant  at  1400  K  may  be  as  high  as  5.63x1  O'10  cm3  s'1. 

The  high  temperature  flowing  afterglow  affords  the  opportunity  to  study  internal 
energy  effects  on  ion-molecule  reaction  kinetics.  Vibrational  excitation  increases  the 
reactivity  of  Ar+  with  02  and  CO  but  hinders  it  for  Ar+  reacting  with  C02.  Only  slight 


vibrational  enhancement  may  arise  for  02+  charge  transfer  to  NO  at  ionospheric 
temperatures.  Obtaining  empirical  values  for  rate  constants  at  high  temperatures  is 
crucial  to  modelers.  Extending  the  technique  to  further  reactions  will  provide  valuable 
insight  into  plasma  chemistry. 
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FIGURES 
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Figure  1:  Rate  constants  for  the  charge  transfer  reactions  of  Ar+  with  02  and  CO 
plotted  against  the  sum  of  the  average  rotational  and  translational  energy.  The  HTFA 
points  are  the  present  results  and  the  FDT  results  are  the  flow  drift  tube  data  of  Dotan 

and  Lindinger.5  The  rate  constants  for  reaction  from  the  v/A>"\  levels  have  been 
calculated  using  the  drift  tube  results  combined  with  the  present  results  assuming  that 
v"=0  and  lof  02  and  CO  react  at  the  drift  tube  rate  and  all  v/JM  react  at  the  same  rate. 
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Figure  2:  Rate  constants  for  the  charge  transfer  reaction  of  Ar+  +  CO2  plotted  against 
the  sum  of  the  average  vibrational,  rotational  and  translational  energy.  The  HTFA 
points  are  the  present  results  and  the  FDT  points  are  the  flow  drift  tube  data  of  Dotan 
and  Lindinger.5 
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Figure  3:  Rate  constants  for  the  charge  transfer  reaction  of  02+  +  NO  plotted  against 
the  temperature  in  K.  The  HTFA  results  are  the  current  results  and  the  FDT  and  FA 
points  are  the  flow  drift  tube  and  flowing  afterglow  results,  respectively,  of  Lindinger  et. 
al.8,9 
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1  Introduction 

This  report  reviews  the  work  done  to  adapt  the  lattice  Boltzmann  (LB)  method 
to  the  prediction  of  atmospheric  events,  including  cloud  development  and  turbulent 
eddies.  These  events  are  relevant  to  battle  site  visibility.  The  scale  on  which  they 
occur  is  in  between  the  macrophysical  and  microphysical  scales  that  are  addressed 
by  the  currently  available  specialized  models.  The  LB  approach  is  selected  because 
it  parallelizes  well  and  therefore  makes  large  computations  feasible. 

The  LB  method  has  proven  useful  for  numerous  fluid  dynamics  problems.  Until 
recently,  the  stable  applications  of  the  method  were  confined  to  isothermal  flows.  Re¬ 
cent  advances  [11]  have  enabled  use  of  the  LB  approach  for  problems  involving  heat 
transfer  without  giving  up  the  stability  of  the  isothermal  version  of  the  method.  It 
should  be  noted  that  while  the  temperature  is  allowed  to  vary,  the  fluid  is  still  de¬ 
scribed  by  an  isothermal  equation  of  state.  This  variable-temperature  improvement  is 
used  as  a  basis  for  additional  modifications  needed  to  address  the  atmospheric  prob¬ 
lems  of  interest.  The  modifications  presented  in  this  report  address  approximating 
an  ideal  gas,  rather  than  an  isothermal  one,  and  simulating  motions  in  a  stratified 
medium.  Experiments  with  parallelization  are  also  reported  here. 

This  report  is  organized  into  sections.  Section  2  briefly  describes  the  lattice  Boltz¬ 
mann  method  and  the  advances  made  by  other  researchers  that  are  relevant  to  the 
atmospheric  problems  of  interest.  The  new  modifications  proposed  by  the  present 
study  are  described  in  Section  3.  Section  4  gives  the  results  and  a  discussion  of  the 
parallelization  experiments. 
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2  The  Lattice  Boltzmann  Method 


A  brief  review  of  the  lattice  Boltzmann  (LB)  method  is  given  here  along  with  a  dis¬ 
cussion  of  the  relevant  improvements  and  applications.  A  more  detailed  presentation 
of  the  basic  method  can  be  found  in  [2]  or  [8].  The  improvements  and  applications 
discussed  here  are  introduced  in  [11],  [12],  and  [13]. 

Perhaps  the  easiest  way  to  understand  the  LB  method  is  to  consider  its  historical 
development.  The  LB  method  was  preceded  by  the  lattice  gas  method  which  uses  in¬ 
tegers  to  represent  particles  that  are  subject  to  translations  and  collisions  to  simulate 
a  flow.  The  particles  occupy  nodes  on  a  regular  lattice  and  are  translated  along  links 
to  neighboring  nodes.  While  this  method  produces  noisy  results  and  some  unphysical 
artifacts,  it  has  the  advantages  of  stability,  parallelizability,  and  a  capability  to  handle 
complicated  boundaries  easily.  The  noise  problem  is  overcome  by  using  real  numbers 
to  represent  particle  distributions  and  computing  collisions  according  to  a  Boltzmann 
equation.  In  particular,  a  distribution  of  particles,  n0(x,  f),  located  at  x  =  (x,  y,  z) 
at  time  t  and  moving  in  direction  a  is  evolved  according  to  the  kinetic  equation 

na(x  +  ea,t-r  1)  =  na(x,t)  +  Qa(x,t),  (1) 

where  ea  is  a  unit  vector  in  direction  a  and  Qa  is  a  collision  operator  for  the  direction 
a  [6].  The  directions  are  aligned  with  the  links  that  connect  neighboring  nodes.  The 
local  (number)  density,  n,  and  the  momentum,  nu  are  recovered  from  the  distribution 
by  the  sums 

n  =  ^na,  nu  =  L^ea.  (2) 

a  a 

This  technique  is  called  the  lattice  Boltzmann  method.  The  LB  collision  operators 
that  were  first  proposed  share  some  of  the  problems  of  the  lattice  gas  method.  Namely, 
these  models  correspond  to  a  momentum  equation  that  contains  an  unphysical  de¬ 
pendence  on  density  in  the  advection  term  and  an  unphysical  dependence  of  pressure 
on  velocity.  These  difficulties  are  resolved  by  the  use  of  the  Bhatnagar,  Gross,  and 
Krook  (BGK)  collision  along  with  a  carefully  defined  equilibrium  distribution  [8]. 

The  BGK  collision  can  be  viewed  as  a  relaxation  method  in  which  the  collision 
produces  na  +  fiQ  =  (1  —  u)na  +  a meaq,  where  u  is  the  relaxation  parameter  and  neaq 
is  the  local  equilibrium  distribution,  which  is  discussed  later.  The  scheme  is  linearly 
stable  for  0  <  u  <  2  with  0  <  u  <  1  resulting  in  under-relaxation  and  1  <  u  <  2 
resulting  in  over-relaxation.  The  corresponding  kinetic  equation  is  usually  written  in 
terms  of  a  relaxation  time,  r  =  as 

na(x  +  ea,  t  +  1)  =  na(x,  t)  -  -[na(x,  t)  -  <9(x,  t)},  (3) 

T 
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and  the  requirement  for  linear  stabiltv  is  r  >  5.  It  will  be  shown  in  a  subsequent 
section  that  the  fluid  viscosity  is  proportional  to  r—  5  so  that  the  stability  requirement 
is  equivalent  to  a  positive  viscosity  requirement. 

There  is  some  flexibility  in  the  choice  of  the  equilibrium  distribution.  For  example 
a  number  of  specifications  for  neq  are  based  on  an  expansion  in  powers  of  the  local 
velocity.  The  coefficients  in  a  truncated  expansion  are  chosen  to  get  the  correct 
behavior  from  the  model.  Typically,  the  goal  is  to  satisfy  the  conservation  of  mass  and 
momentum  equations.  The  unphysical  results  discussed  above  can  be  eliminated  by 
applying  certain  restrictions  to  the  coefficients  in  terms  of  the  number  of  dimensions 
and  the  lattice  structure.  The  details  of  this  process  are  given  in  Section  2.1  for  a 
particular  choice  for  neq. 

It  has  been  demonstrated  that  neq  can  be  determined  so  that  the  LB  method  simu¬ 
lates  continuity,  momentum,  and  energy  equations  [1]  so  that  flow  problems  involving 
heat  transfer  could  be  addressed.  This  approach  incorporates  stability  problems  that 

can  be  avoided  if  only  the  continuity  and  momentum  equations  must  be  satisfied.  An 

» 

alternate  approach  for  addressing  heat  transfer  problems  is  also  available.  It  relies  on 
a  version  of  the  LB  method  that  allows  for  multiple  species,  each  of  which  satisfies 
its  own  continuity  and  momentum  equations  [12].  The  different  species  can  interact 
via  interparticle  forces  and  external  forces.  Heat  transfer  problems  can  be  handled  by 
letting  one  of  the  components  represent  temperature  [11].  This  is  discussed  further 
in  Section  2.2. 

2.1  Equilibrium  Distribution 

This  section  develops  the  requirements  placed  on  the  coefficients  in  an  equilibrium  dis¬ 
tribution  for  a  three-dimensional,  fourteen-direction  model.  To  do  this,  a  Chapman- 
Enskog  expansion  is  used  to  find  the  differential  equations  approximated  by  the  lattice 
Boltzmann  equation  (with  the  BGK  collision).  A  particular  form  for  the  equilibrium 
distribution  is  assumed  and  examined  in  terms  of  the  physics  it  simulates.  The  com¬ 
plete  specification  of  the  equilibrium  distribution  is  given. 

The  lattice  Boltzmann  equation  in  (3)  is  a  discrete  version  of 

na(x  +  6xea,  t  +  5t)  -  na(x,  t)  =  --[na(x,  t)  -  n*9(x.  t)},  (4) 

T 

where  the  local  number  density,  n,  and  velocity,  u,  are  found  from  na  as  shown  in 

(2).  A  Taylor  expansion  of  na  can  be  written 

* 

3  1 

na(x  -r  <Lea,  t  +  5t)  =  rza(x,  t)  +  5xea  •  Vna.+  5t— na  +  -<^eaea  :  VYnana 
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(5) 


-  ,  dna  l,id2 
+Mlea.V—T 

where  everything  on  the  right  hand  side  is  evaluated  at  (x,  t).  The  Chapman-Enskog 
approach  is  used  to  expand  the  time  derivatives  and  the  distribution  function.  Mul¬ 
tiple  time  scales  are  assumed  so  that. 

where  K<  1.  Changes  due  to  convection  are  captured  by  the  faster  (tt)  time  scale, 
while  changes  due  to  diffusion  are  captured  by  the  slower  (f2)  time  scale.  For  near¬ 
equilibrium  conditions,  the  distribution,  na  can  be  expanded  as, 

na  =  Kq  +  +  cn{2)  H - .  (7) 

The  lattice  Boltzmann  equation  given  by  (4)  can  be  rewritten  using  the  Taylor 

and  Chapman-Enskog  expansions,  which  gives  the  leading  order,  0(e),  equation 

> 

e- ' v<"  *  Jinr  =  ■;n*  ’  (8) 

and  the  next  order,  0(e2),  equation 

e  .VfjWj.AnWA/uL  e  ■  VVne?  -  e  .  V-^-r7e<?  -a- r)eq  —  — — (Q) 
a  0  +dh  a  dt2  a  '  2e°ea-  VVn“  '  e°  V5t!  °  '  2dt\  a  -  r  a  •  [  } 

The  leading  order  equation  can  be  used  to  eliminate  the  higher  order  derivatives  so 
that 

<10> 

replaces  the  0(e2)  equation  given  in  (9). 

The  conservation  equations  are  obtained,  to  order  e2,  from  Equations  (8)  and  (10). 
For  example,  by  adding  Equation  (8)  and  the  product  of  e  and  (10),  then  summing 
the  result  over  the  directions,  a,  the  conservation  of  mass  is  obtained  to  order  e2, 


—  -r  V  •  (nu)  =  0.  (11) 

The  momentum  equation  is  found  by  adding  Equation  (8)  and  the  product  of  e  and 
Equation  (10),  multiplying  the  result  by  ea  and  then  summing  over  a,  which  can  be 
rearranged  to  give, 

^H  +  V-  J>aea  nf  +  e  (l  -  =0.  (12) 


Equation  (8)  is  used  to  eliminate  so  that  (12)  is  replaced  by 


— (mi)  4-  V-£>aea 


n”  ~ €(T  -  b  [k + e 


a  V  j  na 


eq 


=  0. 


(13) 


This  is  further  developed  for  a  specific  lattice  and  form  for  flq. 

The  form  adopted  for  the  equilibrium  distribution  is  the  truncated  series 


neaq  =  n(Aa  +  Baea  •  u  4-  Caeaea  :  uu  4-  Dan  ■  u),  (14) 

for  a  three-dimensional  lattice  with  a  =  0, . . . ,  14.  Each  node  is  linked  to  six  of  its 
neighbors,  in  the  north,  south,  east,  west,  upward,  and  downward  directions,  by  links 
of  length  1  and  to  eight  of  its  neighbors  in  the  diagonal  directions  (for  example,  the 
north,  east,  upward  direction)  by  links  of  length  \/3.  Stationary  particles  correspond 
to  a  =  0  with  e0  =  0.  It  can  be  assumed  that  the  coefficients  Aa,  Ba ,  Ca,  and  Da  for 
the  fourteen-direction  model  depend  only  on  |eQ|2  and  the  local  values  for  n  and  u. 
The  requirement,  n  =  X)ana9>  an<3  the  summation  rules  in  the  Appendix  lead  to  the 
coefficient  constraints 


A(0)  +  6A(1)  +  8A(3)  =  1, 
2 C(1)  -r  8C(3)  4-  D{0)  +  6jD(1)  +  8D(3)  =  0, 


(15) 

(16) 


where  the  superscripts  (0),  (1),  and  (3)  refer  to  |ea|2  for  the  given  direction.  In 
addition,  a  distribution  that  satisfies  nu  =  must  have 

Bw  +  4B(3)  =  i  (17) 

jU 

The  remaining  constraints  on  the  coefficients  of  the  equilibrium  distribution  are 
found  by  requiring  that  the  lattice  Boltzmann  equation  is  consistent  with  the  mo¬ 
mentum  conservation  equation.  The  term  V-£aeaeanf ,  from  Equation  (12),  can  be 
simplified  using  the  summation  rules  shown  in  the  Appendix  to  give 


v  •  £  e.e„n? 

a 


where 


2 Vn[A(1)  +  4 A<4>  +  (D{1)  +  4£>(3)  +  4C(3))u  •  u] 

-fl6V  •  (C(3)nuu)  +  2V(C(1)  -  SC^nT,  (18) 


f 

ul 

0 

0 

\ 

T  = 

0 

ul 

0 

{ 

0 

0 
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Uz 
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This  can  be  further  simplified  and  used  to  rewrite  the  0(1)  terms  in  Equation  (12) 
combined  with  the  conservation  of  mass  equation,  (11),  to  obtain  e 

0  =  n-^Al6C(3)nu-Vu4-Vp  +  (16C(3)-l)uV-(mi) 

4-16nu(u  •  VC(3))  4-  V  •  [2 n{C{1)  -  8C(3))T]  4-  0(e),  (20) 

where 

p  =  2n[.4<1)  4-  4A(4>  4-  (4C(3)  +  D{1)  4-  4D(3))u  •  u].  (21) 

A  comparison  of  Equation  (20)  and  the  Euler  equations  leads  to  several  more 
constraints  on  the  coefficients.  In  order  to  eliminate  the  dependence  of  p  on  u,  the 
requirement, 

4 C(3)  4-  £>(1)  4-  4D(3)  =  0,  (22) 

is  imposed,  which  leaves 

p  =  2n(A(1)  4-4A(3)).  (23) 

» 

> 

The  convection  term,  nu  •  Vu,  has  the  required  coefficient  when 

C<3>  =  (24) 

The  last  term  in  Equation  (20)  is  eliminated  by  setting  C M  —  8C^  =  0  so  that 

cm  =  \.  (25) 

With  (22)  -  (25),  Equation  (20)  becomes 

9u 

n— -r  nu  •  Vu  4- Vp  =  O(e).  (26) 

The  order  e  terms  of  Equation  (12)  are  considered  in  order  to  obtain  the  viscous 
terms  in  the  equation  of  motion  and  additional  coefficient  constraints.  Using  the 
form  for  given  by  (14)  and  the  coefficient  constraints  given  above,  these  terms  are 
written 

V  •  Y,  e“ea^  +  e“  '  v)n?  =  JjT  [VP  +  V  '  (nuu)l 
4-16VV  •  CB(3W)  4-  8V  •  V(5(3)nu)  +  2D[(5(1)  -  8B(3))nu]  (27) 

where 

/  d2/dx 2  0  0  \ 

D  =  0  d2/dy 2  0  .  (28) 

VO  0  d2/dz 2 
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The  last  term  in  (27)  contributes  to  the  viscosity  for  certain  one-dimensional  problems 
and  can  be  eliminated  by  setting 

Bw  _  8 B{3)  =  0.  (29) 

Together  with  (17)  this  gives  B ^  =  1  and  B (3)  =  dj.  When  these  coefficients  are 
used  in  (27),  it  can  be  rearranged  using  the  continuity  and  momentum  equations  so 
that  to  order  e 

-  ('  -  5)  V-£  e«e«  (JfT.  +  e„  •  vi  n”  =  -n v!(au)-«( v  nui- i  t  -  i)  Jk(nu) 

(30) 

where  the  coefficients  of  viscosity  and  bulk  viscosity  are 

v  =  \  (r  “  and  Tj  =  7  _  2(A(1)  +  4A(3))  (r  “  .  (31) 

respectively. 

The  coefficients  A^  and  D ^  are  as  yet  undetermined  (as  are  B ^  and  C^\ 
however  these  are  irrelevant  since  they  are  multiplied  by  e0  =  0).  Two  equations 
for  the  three  coefficients  are  provided  by  (16)  and  (22).  Adopting  the  choice 
D used  by  Qian  et  al.  [8]  leads  to  D ^  =  —1  and  D ^  Furthermore, 

the  choice  A(0)  =  |  and  ^  =  |,  also  used  by  Qian  et  al.,  implies  that  A^  =  1  and 
A^3)  =  Note  that  equating  £  to  a  constant  sets  the  sound  speed  to  an  isothermal 
one,  which  in  this  case  has  the  value 

With  these  choices  for  the  equilibrium  distribution  coefficients,  the  equations  of 
motion  in  (11)  and  (13)  can  be  rewritten  as  the  Navier-Stokes  equations 

^-rV-(pu)  =  0,  (32) 

5u  „ 

p—  -r  pn-'Vu  =  -Vp  +  yV2(/?u)  -1-  7?V(V  -  pu),  (33) 

where  p  =  mn  with  m  representing  the  mass  per  particle  and  p  redefined  in  terms  of 

p  (instead  of  n)  so  that 

V  =  |  ■  (34) 

2.2  Passive  Scalar  Temperature 

The  use  of  two  species  of  particles  with  one  representing  mass  and  the  other  repre¬ 
senting  temperature  was  proposed  in  [11]  and  is  briefly  described  here.  The  method 
is  based  on  the  multiple-species  method  proposed  in  [121.  Each  species,  a ,  has  its 


/ 


own  molecular  weight.  ma,  density,  pa  =  m^rSc\  and  relaxation  time,  ra.  The  num¬ 
ber  density  and  velocity  for  each  species,  and  uff,  respectively,  can  be  computed 
independently.  The  two  species  are  assumed  to  have  a  common  velocity, 


Trr 


(35) 


This  common  velocity  may  be  incremented  by  different  amounts  for  each  species  to 
reflect  accelerations  due  to  the  forces  associated  with  each  species.  In  this  case,  where 
species  2  is  temperature,  m2  =  0,  T  =  n the  incremented  common  velocity  is  the 
same  for  both  species.  In  this  way,  the  temperature  moves  with  the  fluid  as  a  passive 
scalar  and  affects  the  fluid  via  a  gravtitational  force.  The  result  is  that,  species  1 
satisfies 


cm 

dt 


dp 

dt 


+  V  •  (pu)  =  0, 


_  Vp  _9 
-r  u  •  Vu  = - h  z/V‘,u  +  g, 


where  p  =  m\rS ^  and  species  2  satisfies 


fdu.VT  =  kV2T, 

at 


where  k  =  |(r2  —  5). 


(36) 

(37) 


(38) 


2.3  A  Note  on  Computer  Implementation 

The  methods  described  above  are  programmed  using  FORTRAN  90  and  computed 
sequentially.  Parallelization  is  discussed  later.  The  code  is  tested  using  several  prob¬ 
lems  with  known  solutions.  For  example,  the  single-species  code  is  checked  by  run¬ 
ning  a  ’viscometer’  problem  and  an  acoustics  problem.  The  viscometer  test  problem, 
proposed  in  [5]  uses  the  decay  of  a  Poiseuille  flow  to  measure  the  viscosity.  This  two- 
dimensional  problem  is  run  in  several  orientations  for  several  values  of  r  to  confirm 
that  the  code  is  running  properly. 

Additional  tests  are  done  using  a  one-dimensional  acoustics  problem.  The  small- 
disturbance  assumption  of  acoustics  is  made  so  that  the  Navier-Stokes  equations  can 
be  linearized  using  the  expansions 


p(x,t)  =  p0(x)  +  epi(x,t)  +  •••, 

(39) 

P(x,i)  =  Po(x)  +  epi(x,  t)  H - , 

(40) 

u(x,  t )  =  u0  (x)  -r  eui  (x,  t)  -T - , 

(41) 
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where  the  subscript  0  refers  to  the  undisturbed  state,  the  subscript  1  is  the  leading 
order  approximation  of  the  disturbance,  and  Kl.  These  expansions,  along  with  the 
assumption  of  isentropic  flow  lead  to  the  following  one-dimensional  approximation  for 


d2pl  po  d2px  d2  dpx  f  ^ 

(42) 

where  v  and  77  are  the  viscosity  and  the  kinematic  viscosity,  respectively.  The  periodic 
solution  has  the  form, 


pi  —  Aeat  cos  ( bt  +  <j>)  sin 


(43) 


where  L  is  the  length  of  the  domain,  <j>  is  a  phase  angle, 
determined  by  the  initial  conditions  and 


<3.11  VI 


k 


are 


constants 


so  that  ^  and  v  +  77  can  be  found  from  (44)  and  (45). 

The  accuracy  of  the  lattice  Boltzmann  solution  to  the  one-dimensional  acoustics 
problem  can  be  assessed  by  applying  the  method  for  given  values  of  A ,  k,  and  L  and 
measuring  the  error  in  ^  (the  square  of  the  sound  speed)  and  v  -f  77.  The  phase 
angle  <j)  is  also  found  in  the  process,  although  it  does  not  play  a  part  in  evaluating 
the  lattice  Boltzmann  solution.  The  values  of  a,  b,  and  <p  for  the  lattice  Boltzmann 
solution  are  taken  to  be  those  that  minimize  the  error  in  pi  defined  by 


£(a,  b,  <f>)  =  UMa,  6,  0;  x{,  t,)  -  p[B{xh  tj)}2,  (46) 

tj  Xi 

where  pi{a,  b,  <f>]  Xi,  tj)  is  calculated  using  (43)  and  p^B(x{,  tj)  is  a  lattice  Boltzmann 
result.  The  minimum  value  for  E  corresponds  to 


(dE  dE_  dE\ 

V  da  ’  db  ’  dp  J  ~  0 


(47) 


The  solution  (a,  6,  <p)  can  be  found  using  a  numerical  root-finding  method. 

This  problem  is  used  to  test  the  lattice  Boltzmann  code  for  lattices  with  20,  40, 
60,  and  80  nodes  in  the  ^-direction  over  a  40  time-step  period.  The  values  k  =  2  and 
r  =  .7  are  used  and  L  is  set  to  the  number  of  nodes  in  the  z-direction.  Equation 
(47)  is  solved  using  a  combination  of  steepest  descents,  to  get  a  rough  approximation, 
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and  Newton's  method,  to  get  a  more  accurate  solution.  The  results  for  a  and  b  are 
used  to  obtain  the  observed  values  for  ^  and  v  -f  r?  to  be  compared  to  theoretical 
values  determined  by  the  lattice  and  the  equilibrium  distribution  as  discussed  in 
Section  2.1.  The  errors  in  the  soundspeed  values  decreased  quadratically  with  lattice 
size  and  ranged  from  .1%  (for  the  80-node  lattice)  to  1.2%  (for  the  20-node  lattice). 
The  smaller-order  effects  due  to  viscosity  are  not  predicted  with  quadratic  accuracy, 
however  the  maximum  error  is  only  1%. 

3  Special  Applications  and  Extensions  of  the  lat¬ 
tice  Boltzmann  Method 

Several  modifications  and  adaptations  of  the  lattice  Boltzmann  method  have  been 
developed  in  working  towards  a  model  suitable  for  predicting  atmospheric  events  that 
involve  moist  air.  The  first  topic  addressed  here,  numerical  boundary  conditions,  is 
relevant  to  a  wide  range  of  problems  in  fluid  dynamics.  Next,  modeling  stratification, 
an  important  feature  of  the  atmosphere,  is  discussed.  Stratification  can  be  achieved 
in  the  lattice  Boltzmann  model  by  applying  the  appropriate  external  forces.  This 
is  explored  below  with  an  example  of  a  gravity  wave  problem.  Finally,  there  is  a 
proposed  modification  to  the  lattice  Boltzmann  method  for  simulating  flow  of  an 
ideal,  rather  than  an  isothermal,  gas. 

3.1  Boundary  Conditions 

Numerical  boundary  conditions  must  be  applied  outside  of  the  domain,  at  the  bound¬ 
ary  nodes,  so  that  during  streaming,  the  appropriate  distributions  are  streamed  into 
the  domain.  This  must  be  done  so  that  the  lattice  Boltzmann  equation  is  satisfied  at 
the  boundary.  If  this  requirement  is  ignored  in  assigning  distributions  at  the  bound¬ 
aries,  then  it  is  effectively  like  using  different  values  of  r  at  these  sites.  The  result 
would  be  a  nonuniform  and  incorrect  viscosity  field.  Dirichlet  and  no-slip  boundary 
conditions  are  addressed  here.  The  method  described  below  is  adapted  from  [11]  for 
general  boundary  shapes. 

For  this  discussion,  it  is  assumed  that  the  initial  condition  assigns  na  =  0  for  any 
node  outside  the  domain  of  the  flow.  Thus,  after  streaming,  nodes  adjacent  to  the 
interior  of  the  domain,  or  boundary  nodes,  have  nonzero  distributions  directed  away 
from  the  domain.  These  directions  are  denoted  a(— ).  The  distributions  directed 
into  the  domain,  or  in  the  a(+)  directions,  are  specified  as  discussed  below.  The 
distributions  that  neither  direct  particles  into  or  out  of  the  domain  have  directions 
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a(0).  These  are  unspecified  and  will  remain  that  way  since  they  are  not  required  for 
computation  purposes. 

The  Dirichlet  condition  in  a  .D-dimensional  domain  specifies  the  number  density,  n, 
and  the  velocity,  u,  at  each  boundary  node,  placing  1  +  D  restrictions  on  the  values 
assigned  to  ria(+)  and  nB(o).  For  example,  consider  a  three-dimensional,  fourteen- 
direction  model  applied  to  a  case  with  a  horizontal  boundary.  The  Dirichlet  condition 
introduces  only  four  equations,  to  evaluate  ten  unknowns.  A  simple  way  of  imposing 
additional  requirements  for  a  general  boundary  shape  follows. 

Conservation  of  mass  and  momentum  are  satisfied  when  the  collision  operator, 

— 1  /r(na  —  neq )  is  applied.  Summing  over  all  directions,  the  collision  results  in 
/  \ 

»?)  +  ZK-  O  +  E(».  -  <’)  =  o.  (48) 

r  W)  «(-)  “(0)  / 

The  a(0)-direction  distributions  can  be  assigned  (somewhat  arbitrarily)  so  that, 

E  n.  =  E  ,  (49) 

a(0)  a(0) 

A  choice  for  nQ(+)  that  satisfies  (48)  is 


in  which  ea(+)  and  e„(_)  have  opposite  directions.  Once  this  assignment  is  made  for 
all  links  directed  into  the  domain,  the  collision  calculation  can  be  carried  out  for  those 
links  using  the  given  density  and  velocity  at  the  boundary.  Note  that  the  approach 
can  be  applied  to  any  boundary  node,  regardless  of  boundary  orientation  or  shape.  An 
alternative  to  selecting  e0(_)  so  that  it  is  opposite  ea(+)  is  to  make  the  selection  based 
on  a  reflection.  The  reflection  rule  would  not  work  for  certain  irregular  boundaries, 
since  it  is  possible  to  reflect  particles  with  an  a(— )  direction  and  obtain  a  new  a(— ) 
direction  rather  than  an  a(+)  direction. 

A  boundary  condition  that  specifies  u  but  not  n  corresponds  to  a  no-slip  boundary. 
This  can  be  implemented  using  a  similar  approach.  A  difference  is  that  in  the  no-slip 
case  some  symmetry  is  assumed  and  the  distributions,  na(+),  are  assigned  using  na(_) 
according  to  a  reflection  or  a  bounce-back  (opposite)  rule.  Therefore,  in  addition  to 
Equations  (48)  and  (49),  the  condition  Tla(-)  na  =  Z!o(+)  na  along  with  u  =  0  gives 


2  Hq(+)  Kg 
Sa(-r)  -da  ~  Sa(— ) 


The  collision  is  then  computed  using  n  as  given  in  Equation  (51)  and  u  =  0  to 
determine  n*^+)  and  using  the  reflected  or  bounced  value  for  na(+)  (as  the  prior-to- 
collision  value). 
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3.2  Gravity  Waves  and  the  Relevant  External  Forces 


A  medium  which  has  as  its  steady  state  a  uniform  density  and  a  potential  temperature 
that  varies  exponentially  with  height  is  discussed  here.  These  conditions  are  relevant 
to  the  problem  of  gravity  waves  in  a  shallow  volume.  The  external  forces  required  to 
maintain  the  desired  stratified  flow  are  presented  below.  Following  that  is  an  outline 
of  a  gravity-wave  problem  along  with  computational  results. 


3.2.1  External  Forces 

External  forces  are  incorporated  into  the  multiple-species  model  to  simulate  the  effects 
of  gravity  based  on  the  approach  in  [13].  This  leads  to  a  general  form  of  the  mass 
balance  equation  for  each  species.  Each  species,  a,  is  acted  upon  by  its  own  external 
force,  F^.  =  n^<T'>ma ga,  where  and  ma  are  the  number  density  and  mass  per 
particle,  respectively  of  species  o  and  gg  is  the  acceleration  due  to  the  force  on 
species  a.  The  model  also  allows  for  interparticle  forces  and  a  distinct  equilibrium 
distribution  function  for  each  species,  however,  these  capabilities  are  not  used  here. 

The  present  application  requires  a  two-species  model,  where  species  1  represents 
the  density,  p  =  m^1),  of  a  fluid  and  species  2  represents  its  potential  temperature, 
0  =  The  molecular  weights  are  mi  =  1,  m2  =  0.  The  external  force  Fx  is 

a  direct  simulation  of  pg,  where  g  is  the  acceleration  due  to  gravity.  The  force  F2 
acts  to  maintain  the  exponential  variation  of  species  2  with  height.  The  forces  are 
specified  in  accordance  with  the  mass  balance  equation  from  [13].  For  species  1  this 
leads  to  the  required  conservation  equation,  regardless  of  the  values  assigned  to  gi 
and  g2.  For  species  2  this  leads  to 

-q£~  +  V  '  (n(2>u)  =  v  +  £>2Vn(2)  +r2n(2)(g!  -  g2)j  ,  (52) 

where  Vg  is  proportional  to  (ra  —  A).  The  assignment  of  gi  is  discussed  first,  then  g2 
is  found  to  satisfy  (52). 

The  density  used  in  the  calculation,  Fx  =  pg,  is  approximated  using  a  Boussinesq 
assumption  as  discussed  in  [11].  It  is  redeveloped  here  for  the  case  where  species 
2  is  the  potential  temperature  (rather  than  the  temperature).  To  establish  how  the 
density  varies  with  the  other  variables,  consider  the  potential  temperature  in  the  form 


9  =  k- 


where  7  is  the  ratio  of  specific  heats  and  k  is  a  constant.  This  is  rewritten  with 
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following  expansions 


6  =  9q(z)  +  #i(z,  y,z,t)  - ,  (54) 

P  =  Po(*)  +Pi(a:,y,2,i)  4 - ,  (55) 

P  =  po  +  pi(x,y,z,t)  H - ,  (56) 

where  0i/#O)  Pi/Po,  and  pi/po  are  at  least  an  order  of  magnitude  smaller  than  one 

(though  not  necessarily  comparable  in  order  to  one  another).  Substituting  the  ex¬ 
pansions  into  (53)  leads  to 

£  =  <5"> 
do  7  Po  Po 

neglecting  higher  order  terms.  Changes-in.  density  are  primarily  due  to  changes  in 
potential  temperature  so  that 

<58> 

Po  t/0 

and  the  density  can  be  written 

p  ^  po  (l  -  tj-^J  .  (59) 

In  the  static,  conductive  state  the  vertical  momentum  equation  becomes 

-Vp  +  pg  =  0.  (60) 

This  can  be  rewritten  using  (59)  as  — Vp*  +  po#i/#og  =  0,  where  p'  has  been  substi¬ 
tuted  for  p+  p0gz  with  g  —  |g|.  Eliminating  the  unknown,  9\,  gives 

Fi  =  ~Po  S  =  Po  9ek  =.  n(1)migi  (61) 

so 

9  -  90 


The  value  of  g2  is  chosen  to  get  the  required  variation  of  species  2  with  height  in 
static  equilibrium.  Species  2  varies  according  to 

n(2)  =  0o(O)  exp  j ,  (63) 

where  &q(0)  and  k  are  constants.  This  state  is  maintained  by  making  the  following 
assignments 

r  =  Ti=  r2, 

V  =  Vx  =  V2, 

Vk 
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(64) 

(65) 

(66) 


where  where  V  =  i(r-  |)  and  e*  is  the  unit  vector  in  the  positive,  vertical  direction. 
Note  that  if  species  2  were  uniform  or  varied  linearly,  (52)  would  reduce  to  the  required 
conservation  of  mass  equation  provided  that  gx  =  g2. 

In  practice  Fx  and  F2  are  applied  to  both  species  because  species  2  is  a  passive 
scalar  and  its  velocity  is  equal  to  the  velocity  of  species  1.  The  forces  are  applied  by 
incrementing  the  (common)  velocity  according  to 


uf-u  +  r 


(67) 


The  equations  of  motion  are 


!  +  V.(pu) 
~  t  u • Vu 

f  +  u-V* 


V p  _ 2 

- T^V2U-rg, 

P 


=  vV2Q. 


(68) 

(69) 

(70) 


These  look  like  Equations  (36)  -  (38),  except  here  the  density  of  species  2,  is  the  po¬ 
tential  temperature,  6,  and  the  thermal  conductivity  is  equal  in  value  to  the  viscosity 
due  to  the  requirement  iq  =  r2. 


3.2.2  A  Gravity- Wave  Problem 

A  gravity-wave  problem  is  done  to  illustrate  the  use  of  the  forces  described  above. 
More  details  on  this  and  related  problems  are  given  in  [4],  [7],  and  [9].  Gravity 
waves,  or  buoyancy  waves,  cause  mountain  lee  waves  and  clear  air  turbulence  and  are 
important  in  atmospheric  modeling. 

For  the  simple  gravity-wave  problem  considered  here,  the  equations  of  motion  can 
be  approximated  so  that  an  exact  solution  is  found.  In  keeping  with  the  Boussi- 
nesq  approximation  used  to  obtain  external  forces,  the  atmosphere  is  assumed  to  be 
incompressible,  so  that  conservation  of  mass  becomes 
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These  equations  are  linearized  for  small  motions  relative  to  the  overall  height  of  the 
atmosphere.  This  corresponds  to  small  changes  in  the  measured  physical  quantities 
and  a  large  Froude  number, 

where  U  is  the  characteristic  velocity,  N  is  the  characteristic  buoyancy  frequency 
(discussed  later),  and  H  is  the  characteristic  vertical  displacement  of  air.  A  two- 
dimensional  problem  is  considered  and  the  physical  quantities  are  expanded  as  in 
(54)  -  (56)  with  the  velocity  u  =  (tt,  w )  expanded  as 

u  —  u0  +  Ui(x,  z,  t)  -i - ,  (75) 

w  =  wi(x,z,t)-i - ,  (76) 

where  the  leading  order  quantities  have  the  subscript  0  and  the  corrections,  which 
are  an  order  of  magnitude  smaller  have  the  subscript  1.  Note  uQ  is  constant.  Using 
the  expansions,  Equation  (72)  gives  the  leading  order  vertical  momentum  equation 
as  the  static  equilibrium  condition 

dpo  f~~7\ 

where  g  =  |g|,  so  that  p0(z )  =  po(0)  —  Pogz-  The  next  order  momentum  equations 
are 

( d  d  \  1  .  _ 


d  d 

m+u°8i 


,  1  dpi 
ui  4 - 

pQ 


{di  +  Uote)Wl  +  70^  ~T09  =  °*  (79) 

where  Equation  (58)  has  been  used  to  replace  p\  in  the  gravity  term.  The  continuity 
equation,  (71),  is  approximated  as 

du-i  dw-i 

^+ir=°  <8°) 

and  the  potential  temperature  equation  (73)  becomes 


dt  Uodx 


6\  +  Wi-p-  =  0. 

CLZ 


These  equations  can  be  rearranged  to  obtain  a  single  equation  with  one  unknown. 
This  is  done  by  taking  a  z  derivative  of  the  x-momentum  equation  and  an  x  derivative 
of  the  z-momentum  equation  and  combining  the  results  to  eliminate  pressure.  This 


f  dux 

dwA 

,  g  od  x 

dx  ) 

\dz  ' 

dx  ) 

90  dx 
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(82) 


The  operator  is  applied  to  this  result.  The  the  term  can  be 

replaced  using  the  continuity  equation  given  by  (80)  and  the  term  4-  u0^)  #1  can 
be  replaced  using  Equation  (81).  The  final  result  is 


dt  :  U°dx 


d^vh^d^wA  ,  Ar2d2Wx 
dx 2  '  dz2  j  "r  dx2 


where 


N2=g 


d(ln#0) 


the  square  of  the  buoyancy  frequency,  is  assumed  to  be  constant.  Equation  (83)  has 
wave  solutions  of  the  form 

w\  =  We1*,  W  =  Wr  +  iWi,  o  —  kx  +  mz  —  vt  (85) 

where  k  and  m  are  wave  numbers  and  v  is  a  frequency.  A  dispersion  relation, 

(i k 2  -I-  mr)(v  —  U(jk )2  —  Ndk2  =  0  (86) 

is  found  by  substituting  the  solution  (85)  into  (83). 

The  streamline  at  the  surface  provides  further  information  about  the  values  of  k 
and  m.  Once  W\  is  found,  Equations  (59)  and  (77)  -  (80)  can  immediately  be  solved 
to  give 

_ 

g(v  -  u0k) 

u  —  uQ  —  ~j^w^  (88) 

s  -  M0)  t1 + jSk) exp  (?)  ’  (89) 

771 

V  —  Po(0)  -  pogz  +  p0—(v  -  u0k)wi.  (90) 

The  real  parts  of  Equations  (85)  and  (87)  -  (90)  correspond  to  the  physical  solutions. 

Consider  steady  gravity  waves  (so  0  =  kx  +  mz)  due  to  a  corrugated  topography 
where  the  elevation  varies  like  cosd.  Provided  that  mz  <  kx  on  the  surface,  the 
shape  of  the  mountain  determines  k  and  using  (86)  the  value  of  m  is  given  by 

N2 

m2  =  —z —  k2.  (91) 


If  m  is  real,  then  there  is  no  decay  with  height.  An  example  of  such  a  solution  is 
shown  in  Figure  1.  The  figure  shows  constant-#  lines,  or  streamlines,  of  the  flow  past 
a  section  of  a  corrugated  surface  containing  a  single  mountain.  The  relevant  constants 
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Figure  1:  Streamlines  for  flow  past  a  corrugated  surface.  The  solid  lines  show  the 
results  of  a  lattice  Boltzmann  computation;  the  dashed  lines  represent  the  exact 
solution  to  the  approximate  gravity-wave  problem. 


are  Jfc  =  §,  m  =  .08,  N2  =  2  x  10"6,  u0  =  .00893,  WT  =  0,  Wt  =  -.000551.  The  solid 
lines  show  the  lattice  Boltzmann  result  after  45,000  time  steps.  The  initial  condition 
is  undisturbed  flow  (with  p  —  po ,  u  =  uq,  w  =  0,  9(z)  =  9o(2))  and  boundary 
conditions  are  Dirichlet  (as  discussed  in  Section  3.1)  at  the  bottom,  Neumann  at  the 
top,  and  periodic  on  the  sides.  The  values  r  =  .52  and  g  =  .003  are  also  used.  The 
dashed  lines  show  the  solution  according  to  Equation  (89).  Causes  of  discrepancies 
include  the  use  of  a  viscous,  nonlinear  model  for  an  inviscid,  linear  problem.  The 
viscosity  in  the  numerical  solution  is  v  =  .00667  and  the  degree  to  which  the  problem 
is  nonlinear  correlates  to  the  Froude  number. 

In  order  to  get  the  Froude  number  for  this  example,  the  characteristic  vertical 
displacement  must  be  written  in  terms  of  the  known  quantities  so  that  Equation  (74) 
can  be  applied.  The  change  in  height  of  the  fluid  moving  along  the  streamline  given 
by  the  mountain  surface  is  /  w{x,  U(x))dt  where  the  surface  is  given  by  z  =  U(x).  For 
a  corrugation  with  period  y,  the  time  required  for  fluid  to  move  from  the  lowest  to 
the  highest  point  is  about  where  u  is  a  characteristic  horizontal  velocity.  Letting 
w  be  the  characteristic  vertical  velocity,  the  characteristic  verticle  displacement  of  air 

H  =  “  (92) 

ku 

For  the  example  of  Figure  1  u  and  w  may  be  taken  to  be  uq  and  \W\,  respectively,  so 
the  Froude  number  is  Fr  —  7.0.  Recall  that  the  linear  solution  is  valid  for  Fr  1. 

Additional  experiments  with  gravity-wave  problems  could  be  done.  With  the 
problem  discussed  above  the  solution  could  be  computed  for  increasingly  finer  grids 
to  check  that  the  ’more  accurate’  solutions  approach  the  solution  to  the  linear  prob¬ 
lem.  Alternatively,  the  solution  could  be  computed  for  increasingly  greater  Froude 
numbers  to  confirm  that  the  agreement  between  the  linear  solution  and  the  LB  so¬ 
lution  improves  as  the  linear  approximation  gets  closer  to  the  full  solution.  A  more 
realistic  lower  boundary  condition  could  be  used  so  that  the  fluid  is  subject  to  a 
no-slip  condition  along  a  mountain  (defined  in  terms  of  lattice  nodes).  Nonlinear 
gravity-wave  problems  can  also  be  solved. 


3.3  A  Modified  Equation  of  State 

The  lattice  Boltzmann  method  simulates  a  fluid  with  an  isothermal  sound  speed. 
While  the  passive  scalar  temperature  method  allows  changes  in  temperature  to  affect 
the  external  forces  on  the  fluid,  the  pressure  remains  a  constant  multiple  of  the  density. 
The  value  of  the  constant  is  determined  by  the  constants  chosen  for  the  coefficients, 
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A(1)  and  A(3)  in  the  equilibrium  distribution.  A  more  general  approach  to  selecting 
A^  and  A^  is  to  let 

/i(x,  t)  =  — .  (93) 

n 

With  =  |  the  remaining  coefficients  are 

-4<" =  \  G  - h)  and  4<3)  =  h  (-5 +  3h)  ■  (94) 

Particle  distributions  should  be  positive.  A  conservative  requirement  that  assures 
this  is  A^,  A*3*  >  0.  This  leads  to  the  condition 

h<h  4  <95> 

To  simulate  ideal-gas  behavior,  h  is  made  proportional  to  the  temperature, 

h  =  RT ,  (96) 

where  the  constant  R  can  be  taken  to  be  the  molecular  weight  multiplied  by  the 
specific  gas  constant. 

Several  test  runs  were  made  to  experiment  with  the  use  of  nonuniform  and  variable 
h.  In  the  two  tests  described  below  small  pressure  jumps  move  through  a  medium. 
Gravity  and  viscosity  can  be  neglected  and  the  acoustic  approximation  is  made,  so 
that  the  disturbances  travel  at  the  soundspeed.  The  situation  can  be  modeled  by  the 
integral  conservation  equations,  simplified  to  reflect  the  assumptions  listed  above. 
The  results  of  this  model  are  compared  to  the  results  of  the  LB  model. 

The  first  test  simulates  two  uniform  media,  such  as  air  and  water,  separated  by  a 
boundary  surface.  Each  medium  has  its  own  density,  px  or  fa.  The  media  are  at  rest 
and  the  pressure,  p0,  is  uniform  so  that  the  boundary  remains  at  rest  and  the  value 
of  h  is  different  in  the  different  media.  An  acoustic  disturbance,  with  pressure  pd  and 
speed  udl  is  initiated  away  from  the  boundary  and  then  passes  through  one  medium 
and  into  the  second.  When  the  disturbance  reaches  the  boundary,  a  reflected  wave  is 
returned  and  a  transmitted  wave  propagates  into  the  second  medium.  The  behavior 
of  the  wave  is  determined  by  the  acoustic  impedance  in  each  medium,  7 =  pjC,, 
where  c,  is  the  soundspeed  in  medium  i.  The  difference  between  the  final  pressure 
and  the  undisturbed  pressure  is  2(pd  -  Po)7?.2/(7L1  -f  7l2)-  The  computations,  using 
the  lattice  Boltzmann  code  for  a  case  in  which  the  densities  of  the  two  media  initially 
differ  by  about  a  factor  of  two,  gives  results  for  density  that  are  within  3%  of  the 
expected  values. 

A  second  test  involves  a  left-traveling  pressure  jump  and  a  right-traveling  pressure 
jump  in  a  single  medium.  The  pressure  jumps  are  small  enough  (less  than  5%)  to  be 
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considered  using  an  acoustic  approximation  and  p/p  is  nearly  constant.  Therefore, 
this  test  problem  could  also  be  done  using  a  lattice  Boltzmann  model  based  on  a  fixed 
value  of  h.  Each  moving  pressure  jump  causes  an  initially  uniform  state  to  jump  to  a 
new  state.  The  two  jumps  collide  and  form  a  new  state  and  two  new  jumps  as  shown 
in  Figure  2.  The  pressure  in  the  new  state  is  within  .3%  of  the  value  given  by  the 
integral  conservation  equations.  While  this  does  not  demonstrate  a  new  capability 
of  the  lattice  Boltzmann  method,  it  does  show  a  successful  approximation  using  the 
variable  h  model. 

It  is  important  to  note  that  the  method  described  does  not  conserve  energy.  A 
simple  example  illustrates  the  difference  in  behavior  of  an  energy-conserving  fluid  and 
a  fluid  that  obeys  the  LB  model  with  variable  h.  If  two  equal  volumes  of  stationary 
fluid  with  differing  densities,  pressures,  and  temperatures  are  allowed  to  mix,  then 
when  equilibrium  is  reached,  the  fluid  is  stationary  and  has  a  new  uniform  state. 
Conservation  of  mass  requires  that  the  new  density  be  the  average  of  the  two  original 
densities.  Conservation  of  energy  for  a  perfect  gas  requires  that  the  new  temperature 
be  a  mass-weighted  average  of  the  original  temperatures.  The  lattice  Boltzmann 
method  with  h  set  according  to  the  perfect  gas  law  will  correctly  predict  the  final 
density,  but  will  incorrectly  predict  the  final  temperature  as  the  average  of  the  original 
temperatures. 

4  Parallelization 

The  efficiency  of  the  lattice  Boltzmann  (LB)  method  is  realized  when  the  computa¬ 
tions  axe  done  using  parallel  processing.  The  opportunity  to  parallelize  the  LB  code 
was  used  to  investigate  a  new  programming  language  called  ZPL.  It  is  a  machine- 
independent,  array  language  developed  at  the  University  of  Washington  in  Seattle, 
Washington.  An  advantage  of  ZPL  is  that  a  single  program  can  be  compiled  for 
different  single-  or  multiple-processor  machines  by  using  a  machine-specific  compiler 
designed  to  exploit  the  target  computer.  Upon  compilation  the  ZPL  code  produces 
C  code,  which  is  then  compiled  by  the  local  C  compiler.  When  a  ZPL  program  is 
compiled  for  sequential  computations  the  performance  of  the  code  is  comparable  to 
that  of  a  C  program  written  for  a  single  processor.  A  drawback  to  ZPL  is  that  it  is 
new  and  still  under  development  so  that  certain  features  are  not  yet  available  and 
presumably  the  chances  of  discovering  a  bug  are  greater  for  ZPL  than  for  a  more 
established  compiler. 

To  experiment  with  ZPL  the  basic  LB  method  is  coded  using  ZPL  and  the  pro- 
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Figure  2:  '.Moving  pressure  jumps.  The  solid  curve  shows  two  pressure  jumps  (traveling 
towards  the  center)  before  they  collide,  at  time=30,  the  dashed  curve  shows  the  new 
pressure  jumps  at  time=50  (traveling  away  from  the  center)  produced  by  the  collision. 


gram  is  compiled  and  run  on  different  machines.  The  program  solves  the  viscometer 
problem  on  a  10  x  10  x  10  grid  over  500  time  steps  (so  that  a  total  of  5  x  10°  indi¬ 
vidual  LB  computations  are  done).  The  first  machine  used  is  an  SGI  Origin  200,  a 
two-processor,  180  megahertz  machine.  When  the  code  is  compiled  and  run  sequen¬ 
tially  the  run  time  required  is  1:56  (given  as  minutes:seconds).  The  run  time  for  the 
same  computation  done  in  parallel  is  1:01.  A  more  detailed  look  at  the  run  times 
shows  that  the  speed-up  is  primarily  due  to  parallelizing  the  collision  computations 
which  approximately  halves  the  collision  computation  time  from  1:54  to  0:58.  The 
time  used  for  streaming  is  nearly  the  same  for  both  approaches  (the  times  are  0:03  for 
sequential  and  0:04  for  parallel  computations).  This  is  as  expected  since  it  is  the  col¬ 
lision  calculations  that  require  only  local  information  and  therefore  lend  themselves 
to  parallelization.  The  results  also  demonstrate  that  most  of  the  compute  time  is 
accounted  for  by  the  collisions. 

The  second  machine  used  to  try  the  ZPL  program  is  an  SGI  Onyx,  a  four- 
processor,  200  megahertz  machine.  The  very  same  code  that  worked  on  the  Origin 
gave  incorrect  results  on  the  Onyx.  Experiments  with  several  simpler  test  codes  and 
many  exchanges  with  the  ZPL  support  group  at  the  University  of  Washington  point 
to  a  bug  in  the  ZPL  compiler  for  the  Onyx.  While  the  specific  problem  has  not  been 
identified,  the  ZPL  support  group  is  willing  to  work  on  it  given  a  short  program 
that  isolates  the  problem.  Producing  such  a  code  and  resolving  the  problem  would 
be  called  for  if  ZPL  use  is  to  be  continued.  This  matter  is  being  considered.  In 
summary,  when  ZPL  works  correctly,  it  performs  as  advertised,  but  its  reliability  is 
questionable. 

5  Conclusion 

The  modifications  proposed  by  the  present  study  are  elements  of  an  atmospheric 
event  simulator.  In  addition  to  the  extensions  suggested  in  the  previous  sections, 
additional  modifications  are  called  for  in  order  to  simulate  cloud-forming  processes. 
In  particlular  vapor  and  liquid  water  can  be  incorporated  as  two  additional  model 
components.  As  flow  conditions  change,  mass  may  be  transferred  from  one  compo¬ 
nent  to  the  other  to  reflect  phase  changes.  Phase  transitions  can  be  specified  using 
the  cloud  parameterization  method  used  in  the  Penn  State/NCAR  Mesoscale  Model 
(MM5)  [3]. 
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6  Appendix 


Useful  sums  from  [10]  for  evaluating  terms  in  LBE  model  with  on  a  fourteen-direction 


lattice,  where  the  links  have  lengths  1  and  \/3, 

42)  =  (97) 

a 

=  2[u;(l)  +  4w(3)]Sij  (98) 

Eijh  =  5Z^(leal  2)eaieajeakeal  (99) 

a 

=  8u;(3)A^  -I-  2(uf(l)  —  8u/(3)]<^4)  (100) 

where 

A^  =  SijSki  +  SikSji  +  5u5jk,  (101) 

5(4)  =  +5ijSkiSik,  (102) 


and  the  subscripts  (i,  j,  k ,  l )  refer  to  the  coordinate  directions  and  w  is  a  function  of 
the  length  of  link  a. 
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